#### Setting environment ####
require(stringi)
require(MCMCpack)
require(doParallel)
require(mgcv)

# Function to compute posterior means, standard deviations, and credible intervals if specified
posterior.summary <- function(x, CI = TRUE) {
  posterior.mean <- mean(x)
  if (CI == TRUE) {
    HPD.interval <- HPDinterval(mcmc(x))
    c(posterior.mean, HPD.interval)
  } else {
    posterior.sd <- sd(x)
    c(posterior.mean, posterior.sd)
  }
}

#### Data preparation ####
## Voting records
voting.records <- read.csv("voting_records.csv")
voting.records$date <- as.Date(voting.records$date)

# Convert data into a matrix format for subsequent analyses
data.matrix <- as.matrix(voting.records[, -c(1:3)])
rownames(data.matrix) <- 1:nrow(data.matrix)

# Information on nonunanimous cases
nonunanimous.cases <- voting.records[, 1:3]
nonunanimous.cases$year <- as.numeric(format(nonunanimous.cases$date, "%Y"))
nonunanimous.cases$n.justices <- rowSums(! is.na(data.matrix))

## Justice information
justice.data <- read.csv("justice_data.csv")
justice.data$begin <- as.Date(justice.data$begin)
justice.data$end <- as.Date(justice.data$end)
justice.data$end[is.na(justice.data$end)] <- as.Date("2024-01-01")
justice.data$pch <- ifelse(justice.data$career == 1, 21, 
                           ifelse(justice.data$career == 2, 22, 
                                  ifelse(justice.data$career == 3, 4, 
                                         ifelse(justice.data$career == 4, 24, 23))))
justice.data$bg <- c("black", "white", NA, "darkgray", "lightgray")[justice.data$career]

# Calculate the average tenure length of justices in years
as.numeric(mean(justice.data$end - justice.data$begin)) / 365

# List of Chief Justices
chief.justices <- justice.data$name[justice.data$chief == 1]

## Case information
case.data <- read.csv("case_data.csv")

# Combine the voting records and case information
nonunanimous.case.data <- merge(nonunanimous.cases, case.data[, paste0("V", 1:3)], 
                                by.x = "filename", by.y = "V1", 
                                all.x = TRUE, all.y = FALSE, sort = FALSE)

# Classification of civil, criminal, and administrative cases
nonunanimous.case.data$type <- ifelse(stri_detect_regex(nonunanimous.case.data$V2, "[行|分]"), "administrative", 
                                      ifelse(stri_detect_regex(nonunanimous.case.data$V2, "[オ|受|許|ク|テ|マ|ヤ]"), "civil", "criminal"))

## Number of observations
# Number of cases (including multifaceted cases) and justices
dim(data.matrix)

# Numbers of Grand Bench and petty bench cases
table(nonunanimous.case.data$n.justices > 5)

# Number of non-missing observations
sum(! is.na(data.matrix))

# Number of unique cases
length(unique(nonunanimous.case.data$filename))

# Number of multifaceted cases
table(table(nonunanimous.case.data$filename))

# Number of cases in each year (Figure A1)
year.table <- table(factor(nonunanimous.cases$year, 1948:2023))
year.table.GB <- table(factor(nonunanimous.cases$year[nonunanimous.cases$n.justices > 5], 1948:2023))

tiff("Figure_A1.tiff",
     width = 5, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(4, 4, 1, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(1948, 2023), ylim = c(0, 35), 
     xlab = "Year", ylab = "Number of cases", xaxt = "n", yaxt = "n")
abline(v = seq(1948, 2023, 10), col = "gray")
abline(h = seq(0, 35, 5), col = "gray")
box(bty = "L")
lines(1948:2023, year.table, type = "b", pch = 19)
lines(1948:2023, year.table.GB, type = "b", pch = 21, bg = "white")
legend("topright", legend = c("Total", "Grand Bench cases"), 
       pch = c(19, 21), bg = "white", pt.bg = c(NA, "white"))
axis(1, at = seq(1948, 2023, 10), labels = c(seq(48, 98, 10), "08", "18"), lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

#### Preliminary estimation of the static IRT model (Appendix B.2) ####
## Estimate the static IRT model for every 7-8 years period
start.year <- c(seq(1948, 1976, 7), seq(1984, 2016, 8))
end.year <- c(start.year[2:10] - 1, 2023)

cluster <- makeCluster(10)
registerDoParallel(cluster)
static.results <- foreach(i = 1:10, .packages = "MCMCpack") %dopar% {
  data.matrix.sub <- data.matrix[nonunanimous.case.data$year >= start.year[i] &
                                   nonunanimous.case.data$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  data.matrix.sub <- data.matrix.sub[, target.justices.sub]
  MCMCirt1d(datamatrix = t(data.matrix.sub), 
            burnin = 100000, mcmc = 100000, thin = 50, 
            seed = list(NA, i), AB0 = matrix(c(0.64, 0, 0, 0.64), 2, 2))
}
stopCluster(cluster)

# save(static.results, file = "static_results.Rdata")
# load("static_results.Rdata")

## Estimated latent trait parameters for each period (Figure A2)
tiff("Figure_A2.tiff",
     width = 5.6, height = 7, unit = "in", pointsize = 6, res = 1200, compression = "lzw")
layout(matrix(1:10, 5, 2, byrow = TRUE))
par(mar = c(2.5, 0.5, 2, 0.5), lwd = 0.25)
for (i in 1:10) {
  data.matrix.sub <- data.matrix[nonunanimous.case.data$year >= start.year[i] &
                                   nonunanimous.case.data$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  justice.data.sub <- justice.data[target.justices.sub, ]
  results.theta <- static.results[[i]][, stri_detect_regex(colnames(static.results[[i]]), "theta")]
  iteration.mean <- rowMeans(results.theta)
  iteration.sd <- apply(results.theta, 1, sd)
  # Post-process latent trait parameters (see Appendix B.1)
  for (j in 1:nrow(results.theta)) {
    results.theta[j, ] <- (results.theta[j, ] - iteration.mean[j]) / iteration.sd[j]
  }
  posterior.mean <- colMeans(results.theta)
  posterior.CI <- HPDinterval(mcmc(results.theta), 0.95)
  posterior.order <- order(posterior.mean)
  plot(NULL, NULL, bty = "n", xlim = c(-4.5, 4.5), ylim = c(1, length(posterior.mean)),
       main = paste0(start.year[i], "\u2013", end.year[i]),
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = -4:4, col = "gray")
  segments(posterior.CI[posterior.order, 1], c(1:nrow(posterior.CI)),
           posterior.CI[posterior.order, 2], c(1:nrow(posterior.CI)))
  points(posterior.mean[posterior.order], 1:length(posterior.mean),
         pch = justice.data.sub$pch[posterior.order],
         bg = justice.data.sub$bg[posterior.order], cex = 0.6)
  text(ifelse(rep(c(0, 1), 100)[1:length(posterior.mean)] == 0,
              posterior.CI[posterior.order, 1],
              posterior.CI[posterior.order, 2]), 1:length(posterior.mean),
       justice.data.sub$english[posterior.order],
       pos = rep(c(2, 4), 100)[1:length(posterior.mean)], cex = 0.6,
       font = ifelse(justice.data.sub$name[posterior.order] %in% chief.justices, 2, 1))
  axis(1, at = -4:4, lwd = 0.25)
}
dev.off()

## Identify and list the most extreme justices based on their latent trait estimates for each period
extreme.judges.jp <- extreme.judges.en <- list()
for (i in 1:10) {
  data.matrix.sub <- data.matrix[nonunanimous.case.data$year >= start.year[i] &
                                   nonunanimous.case.data$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  justice.data.sub <- justice.data[target.justices.sub, ]
  results.theta <- static.results[[i]][, stri_detect_regex(colnames(static.results[[i]]), "theta.")]
  posterior.mean <- colMeans(results.theta)
  posterior.order <- order(posterior.mean)
  extreme.judges.jp[[i]] <- as.character(justice.data.sub$name[posterior.order][c(1, length(posterior.order))])
  extreme.judges.en[[i]] <- as.character(justice.data.sub$english[posterior.order][c(1, length(posterior.order))])
}
extreme.judges.en

# Signs of the latent traits are fixed in the estimation of the dynamic IRT model
# Figure A2 tells us which side in a period is connected to which side in the next period.
sign <- c(-1, 1, 1, 1, 1, -1, -1, 1, 1, -1)
theta.constraints <- list()
for (i in 1:10) {
  if (sign[i] == 1) {
    theta.constraints[[2 * i - 1]] <- "-"
    theta.constraints[[2 * i]] <- "+"
  } else {
    theta.constraints[[2 * i - 1]] <- "+"
    theta.constraints[[2 * i]] <- "-"
  }
  names(theta.constraints)[2 * i - 1] <- extreme.judges.jp[[i]][1]
  names(theta.constraints)[2 * i] <- extreme.judges.jp[[i]][2]
}

#### Estimation of the dynamic IRT model ####
# Year id
year <- as.numeric(nonunanimous.case.data$year) - 
  as.numeric(min(nonunanimous.case.data$year)) + 1
# Number of justices
N <- ncol(data.matrix)
# Number of cases
J <- nrow(data.matrix)
# Number of years
t <- max(year)

cluster <- makeCluster(10)
registerDoParallel(cluster)
dynamic.list <- foreach(i = 1:3, .packages = "MCMCpack") %dopar% {
  MCMCdynamicIRT1d(datamatrix = t(data.matrix), 
                   item.time.map = year, 
                   theta.constraints = theta.constraints, 
                   burnin = 200000, mcmc = 200000, 
                   thin = 100, seed = list(NA, i), 
                   tau2.start = 0.1, A0 = 0.64, B0 = 0.64)
}
stopCluster(cluster)

# save(dynamic.list, file = "dynamic_list.Rdata")
# load("dynamic_list.Rdata")

# Convergence check
dynamic.results <- mcmc.list(dynamic.list)
which(gelman.diag(dynamic.results, multivariate = FALSE)$psrf[, 1] > 1.1)

## Post-processing of the posterior draws (see Appendix B.1)
raw.MCMC.draws <- normalized.MCMC.draws <- 
  rbind(dynamic.results[[1]], dynamic.results[[2]], dynamic.results[[3]])

# First-year latent trait parameter for each justice
initial.period.column <- rep(NA, N)
for (i in 1:N) {
  initial.period.column[i] <- which(stri_detect_regex(colnames(raw.MCMC.draws), 
                                                      paste0("theta.", justice.data$name[i])))[1]
}

# Convert the raw sampled values into the normalized values
iteration.mean <- rowMeans(raw.MCMC.draws[, initial.period.column])
iteration.sd <- apply(raw.MCMC.draws[, initial.period.column], 1, sd)
for (i in 1:nrow(raw.MCMC.draws)) {
  normalized.MCMC.draws[i, stri_detect_regex(colnames(normalized.MCMC.draws), "theta")] <- 
    (raw.MCMC.draws[i, stri_detect_regex(colnames(raw.MCMC.draws), "theta")] - iteration.mean[i]) / iteration.sd[i]
  normalized.MCMC.draws[i, stri_detect_regex(colnames(normalized.MCMC.draws), "alpha")] <- 
    raw.MCMC.draws[i, stri_detect_regex(colnames(raw.MCMC.draws), "alpha")] - 
    iteration.mean[i] * raw.MCMC.draws[i , stri_detect_regex(colnames(raw.MCMC.draws), "beta")]
  normalized.MCMC.draws[i , stri_detect_regex(colnames(normalized.MCMC.draws), "beta")] <- 
    iteration.sd[i] * raw.MCMC.draws[i , stri_detect_regex(colnames(raw.MCMC.draws), "beta")]
}
normalized.MCMC.draws <- mcmc(normalized.MCMC.draws)

#### Latent trait parameters (justices' ideal points) ####
# Record each case's year
data.year <- t(data.matrix)
for (i in 1:J) {
  data.year[, i] <- ifelse(is.na(t(data.matrix)[, i]), NA, year[i])
}

# Record each justice's first and last years
# These years do not necessarily correspond to justices' tenure (see Footnote 19)
justice.year <- matrix(NA, N, 2)
for (i in 1:N) {
  justice.year[i, 1] <- min(data.year[i, ], na.rm = TRUE)
  justice.year[i, 2] <- max(data.year[i, ], na.rm = TRUE)
}

## Posterior draws for each justice's latent trait parameter
# ideal.point.sample.by.year: posterior draws for each justice-year
# ideal.point.year: posterior mean and credible interval for ideal.point.sample.by.year
# ideal.point.whole: posterior draws for each justice's average over their tenure
# posterior.mean: posterior means for ideal.point.whole
# posterior.CI: credible intervals for ideal.point.whole
# posterior.order: ascending order of posterior.mean
ideal.point.sample.by.year <- array(NA, c(N, t, nrow(normalized.MCMC.draws)))
ideal.point.year <- array(NA, c(N, t, 3))
ideal.point.whole <- matrix(NA, nrow(normalized.MCMC.draws), N)
for (i in 1:N) {
  temp.draws <- normalized.MCMC.draws[, stri_detect_regex(colnames(normalized.MCMC.draws), paste0("theta.", justice.data$name[i]))]
  ideal.point.sample.by.year[i, justice.year[i, 1]:justice.year[i, 2], ] <- temp.draws
  if (is.null(dim(temp.draws))) {
    ideal.point.year[i, justice.year[i, 1]:justice.year[i, 2], 1] <- mean(temp.draws)
    ideal.point.whole[, i] <- temp.draws
  } else {
    ideal.point.year[i, justice.year[i, 1]:justice.year[i, 2], 1] <- colMeans(temp.draws)
    ideal.point.whole[, i] <- rowMeans(temp.draws)
  }
  ideal.point.year[i, justice.year[i, 1]:justice.year[i, 2], 2:3] <- HPDinterval(mcmc(temp.draws))
}
rownames(ideal.point.year) <- justice.data$english
posterior.mean <- colMeans(ideal.point.whole)
names(posterior.mean) <- justice.data$name
posterior.CI <- HPDinterval(mcmc(ideal.point.whole))
posterior.order <- order(posterior.mean)

## Estimated latent trait parameters of justices averaged over their tenures (Figure 1)
tiff("Figure_1.tiff",
     width = 5.6, height = 7, unit = "in", pointsize = 8, res = 1200, compression = "lzw")
par(mar = c(2, 0.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-4.5, 4.5), ylim = c(1, N), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = -4:4, col = "gray")
segments(posterior.CI[posterior.order, 1], c(1:nrow(posterior.CI)), 
         posterior.CI[posterior.order, 2], c(1:nrow(posterior.CI)))
points(posterior.mean[posterior.order], 1:length(posterior.mean), 
       pch = justice.data$pch[posterior.order], 
       bg = justice.data$bg[posterior.order], cex = 0.6)
text(ifelse(rep(c(0, 1), 100)[1:length(posterior.mean)] == 0, 
            posterior.CI[posterior.order, 1], 
            posterior.CI[posterior.order, 2]), 1:length(posterior.mean), 
     justice.data$english[posterior.order], 
     pos = rep(c(2, 4), 100)[1:length(posterior.mean)], cex = 0.6, 
     font = ifelse(justice.data$name[posterior.order] %in% chief.justices, 2, 1))
legend("bottomright", pch = c(16, 22, 4, 24, 23), 
       pt.bg = c(NA, "white", NA, "darkgray", "lightgray"), 
       legend = c("Career Judge", "Attorney", "Prosecutor", "Bureaucrat", "Law Professor"))
axis(1, at = -4:4, lwd = 0.5)
dev.off()

## Estimated year-by-year latent trait parameters of justices (Figure A3)
tiff("Figure_A3_1.tiff",
     width = 6, height = 8, unit = "in", pointsize = 5, res = 1200, compression = "lzw")
layout(matrix(1:48, 8, 6, byrow = TRUE))
par(mar = c(1, 1, 2.5, 0.5), lwd = 0.25)
for (i in 1:48) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(1, t), ylim = c(-5, 5), 
       main = justice.data$english[i], 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = seq(1, t, 10), col = "gray")
  abline(h = seq(-4, 4, 2), col = "gray")
  box()
  segments(1:t, ideal.point.year[i, , 2], 1:t, ideal.point.year[i, , 3])
  points(1:t, ideal.point.year[i, , 1], 
         pch = justice.data$pch[i],  
         bg = justice.data$bg[i], cex = 0.6)
  mtext(c(seq(48, 98, 10), "08", "18"), side = 1, line = -1, 
        at = seq(1, t, 10), cex = 0.6)
  mtext(seq(-4, 4, 2), side = 2, line = -1, at = seq(-4, 4, 2), cex = 0.6)
}
dev.off()
tiff("Figure_A3_2.tiff",
     width = 6, height = 8, unit = "in", pointsize = 5, res = 1200, compression = "lzw")
layout(matrix(1:48, 8, 6, byrow = TRUE))
par(mar = c(1, 1, 2.5, 0.5), lwd = 0.25)
for (i in 49:96) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(1, t), ylim = c(-5, 5), 
       main = justice.data$english[i], 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = seq(1, t, 10), col = "gray")
  abline(h = seq(-4, 4, 2), col = "gray")
  box()
  segments(1:t, ideal.point.year[i, , 2], 1:t, ideal.point.year[i, , 3])
  points(1:t, ideal.point.year[i, , 1], 
         pch = justice.data$pch[i], 
         bg = justice.data$bg[i], cex = 0.6)
  mtext(c(seq(48, 98, 10), "08", "18"), side = 1, line = -1, 
        at = seq(1, t, 10), cex = 0.6)
  mtext(seq(-4, 4, 2), side = 2, line = -1, at = seq(-4, 4, 2), cex = 0.6)
}
dev.off()
tiff("Figure_A3_3.tiff",
     width = 6, height = 8, unit = "in", pointsize = 5, res = 1200, compression = "lzw")
layout(matrix(1:48, 8, 6, byrow = TRUE))
par(mar = c(1, 1, 2.5, 0.5), lwd = 0.25)
for (i in 97:144) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(1, t), ylim = c(-5, 5), 
       main = justice.data$english[i], 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = seq(1, t, 10), col = "gray")
  abline(h = seq(-4, 4, 2), col = "gray")
  box()
  segments(1:t, ideal.point.year[i, , 2], 1:t, ideal.point.year[i, , 3])
  points(1:t, ideal.point.year[i, , 1], 
         pch = justice.data$pch[i], 
         bg = justice.data$bg[i], cex = 0.6)
  mtext(c(seq(48, 98, 10), "08", "18"), side = 1, line = -1, 
        at = seq(1, t, 10), cex = 0.6)
  mtext(seq(-4, 4, 2), side = 2, line = -1, at = seq(-4, 4, 2), cex = 0.6)
}
dev.off()
tiff("Figure_A3_4.tiff",
     width = 6, height = 8, unit = "in", pointsize = 5, res = 1200, compression = "lzw")
layout(matrix(1:48, 8, 6, byrow = TRUE))
par(mar = c(1, 1, 2.5, 0.5), lwd = 0.25)
for (i in 145:N) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(1, t), ylim = c(-5, 5), 
       main = justice.data$english[i], 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = seq(1, t, 10), col = "gray")
  abline(h = seq(-4, 4, 2), col = "gray")
  box()
  segments(1:t, ideal.point.year[i, , 2], 1:t, ideal.point.year[i, , 3])
  points(1:t, ideal.point.year[i, , 1], 
         pch = justice.data$pch[i], 
         bg = justice.data$bg[i], cex = 0.6)
  mtext(c(seq(48, 98, 10), "08", "18"), side = 1, line = -1, 
        at = seq(1, t, 10), cex = 0.6)
  mtext(seq(-4, 4, 2), side = 2, line = -1, at = seq(-4, 4, 2), cex = 0.6)
}
dev.off()

#### Discrimination parameters ####
# Extract draws for discrimination parameters from MCMC results
results.discrimination <- normalized.MCMC.draws[, stri_detect_regex(colnames(normalized.MCMC.draws), "beta.")]

# Calculate the posterior means of the discrimination parameters
discrimination.mean <- colMeans(results.discrimination)

# Compute discrimination power as the absolute values of the discrimination parameters
discrimination.power <- abs(discrimination.mean)

## Identify Grand Bench cases with high discrimination power for inclusion in Table 1
# Note: Only the case with the highest discrimination power is listed in Table 1, among those with similar primary dispute points.
# The primary point of dispute in each case cannot be identified without reference to the judgment texts.
# The case numbers of those displayed in Table 1 are as follows:
# 50906, 51311, 51255, 87091, 56279, 51800, 90412, 56260, 50269, 51371, 50499, 52690, 53752, 51645, 56516, 70667, and 92446
divisive.cases <- data.frame(filename = nonunanimous.case.data[, "filename"], 
                             date = nonunanimous.case.data[, "date"], 
                             case = nonunanimous.case.data[, "V3"], 
                             beta = round(discrimination.mean, 2), 
                             n.justices = nonunanimous.case.data[, "n.justices"])[order(discrimination.power, decreasing = TRUE), ]
rownames(divisive.cases) <- NULL
head(subset(divisive.cases, n.justices > 5), 50)

## Calculate the probability of justices being completely divided on the scale by chance (Footnote 26)
# Five-justice petty bench
round(2 * (5 - 1) / (2 ^ 5 - 2), 3)
# Four-justice petty bench
round(2 * (4 - 1) / (2 ^ 4 - 2), 3)
# Fifteen-justice Grand Bench
round(2 * (15 - 1) / (2 ^ 15 - 2), 3)
# Fourteen-justice Grand Bench
round(2 * (14 - 1) / (2 ^ 14 - 2), 3)

#### Concordance with the anecdote of Ishida court ####
# Extract the ideal points for justices in the Yokota court in their last year of service
# and those of justices in Ishida courts in their first year of service
Yokota.justices <- c("Matsuda, J", "Kusaka, A", "Kido, Y", "Irie, T", "Osabe, K", 
                     "Iimura, Y", "Matsumoto, M", "Iwata, M", "Shimomura, S", 
                     "Irokawa, K", "Tanaka, J")
Ishida.justices <- c("Fujibayashi, M", "Okahara, M", "Ogawa, N", "Shimoda, T", 
                     "Kishi, S", "Amano, B", "Sakamoto, Y", "Kishigami, Y", 
                     "Erikuchi, K", "Otsuka, K", "Takatsuji, M")
Yokota.end <- Ishida.start <- rep(NA, 11)
for (i in 1:11) {
  Yokota.end[i] <- tail(na.omit(ideal.point.year[Yokota.justices[i], , 1]), 1)
  Ishida.start[i] <- head(na.omit(ideal.point.year[Ishida.justices[i], , 1]), 1)
}

## Replacement of justices in the Ishida Court (Figure 2(a))
# Label the figure with classifications from Yamamoto (1994)
Yokota.LC <- justice.data$Yamamoto[match(Yokota.justices, justice.data$english)]
Ishida.LC <- justice.data$Yamamoto[match(Ishida.justices, justice.data$english)]
Yokota.labels <- paste0(Yokota.justices, " (", Yokota.LC, ")")
Ishida.labels <- paste0(Ishida.justices, " (", Ishida.LC, ")")

# Set label positions to avoid overlapping in the figure
Yokota.adjust <- c(0, 0, 0.13, 0, 0, -0.08, 0, 0.07, 0, -0.07, 0)
Ishida.adjust <- c(-0.02, 0, 0, 0.02, 0, 0, 0, 0.13, 0.02, 0, 0.1)

tiff("Figure_2_a.tiff",
     width = 4, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(0, 0, 0, 0), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(0, 3), ylim = c(-2, 2.2), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(1, -2, 1, 2)
segments(2, -2, 2, 2)
points(rep(1, 11), Yokota.end, pch = 19)
points(rep(2, 11), Ishida.start, pch = 19)
segments(1, Yokota.end, 0.75, Yokota.end + Yokota.adjust)
segments(2, Ishida.start, 2.25, Ishida.start + Ishida.adjust)
text(0.75, Yokota.end + Yokota.adjust, Yokota.labels, pos = 2)
text(2.25, Ishida.start + Ishida.adjust, Ishida.labels, pos = 4)
for (i in 1:11) {
  arrows(1.05, Yokota.end[i], 1.95, Ishida.start[i], length = 0.03)
}
text(0.45, 2.2, "M. Yokota Court", font = 2)
text(2.55, 2.2, "Ishida Court", font = 2)
text(1:2, -2, "\u2212", pos = c(2, 4))
text(1:2, 2, "+", pos = c(2, 4))
dev.off()

## Justices’ position in Tokyo-to Kyoso and Zennorin Keishokuho cases (Figure 2(b))
# Row numbers of Tokyo-to Kyoso case (50812) and Zennorin Keishokuho case (50906) in nonunanimous.case.data
Tokyoso.pos <- which(nonunanimous.case.data$filename == 50812)
Zennorin.pos <- which(nonunanimous.case.data$filename == 50906)

# List of justices in these cases
data.matrix.en <- data.matrix
colnames(data.matrix.en) <- justice.data$english
Tokyoso.justices <- names(na.omit(unlist(data.matrix.en[Tokyoso.pos, ])))
Zennorin.justices <- names(na.omit(unlist(data.matrix.en[Zennorin.pos, ])))

# Ideal points of the justices
Tokyoso.ideal.point <- ideal.point.year[Tokyoso.justices, nonunanimous.case.data$year[Tokyoso.pos] - 1947, 1]
Zennorin.ideal.point <- ideal.point.year[Zennorin.justices, nonunanimous.case.data$year[Zennorin.pos] - 1947, 1]

# Cut point in these cases (see Footnote 28)
Tokyoso.alpha <- normalized.MCMC.draws[, paste0("alpha.", Tokyoso.pos)]
Tokyoso.beta <- normalized.MCMC.draws[, paste0("beta.", Tokyoso.pos)]
Tokyoso.cutpoint <- -1 * Tokyoso.alpha / Tokyoso.beta

Zennorin.alpha <- normalized.MCMC.draws[, paste0("alpha.", Zennorin.pos)]
Zennorin.beta <- normalized.MCMC.draws[, paste0("beta.", Zennorin.pos)]
Zennorin.cutpoint <- -1 * Zennorin.alpha / Zennorin.beta

# Figure labels with classification of each justice's opinion
Tokyoso.labels <- paste0(Tokyoso.justices, " (", 
                         justice.data$Tokyoso[match(Tokyoso.justices, justice.data$english)], ")")
Zennorin.labels <- paste0(Zennorin.justices, " (", 
                          justice.data$Zennorin[match(Zennorin.justices, justice.data$english)], ")")

# Set label positions to avoid overlapping in the figure
Tokyoso.adjust <- c(0.05, -0.05, 0.05, 0, -0.05, -0.05, 0.07, 
                    0, -0.07, 0.05, -0.12, -0.01, 0.03, -0.05)
Zennorin.adjust <- c(-0.04, -0.03, 0.06, 0.12, -0.06, 0, -0.12, 
                     0.03, 0, 0, 0, 0, 0.05, 0, -0.04)

tiff("Figure_2_b.tiff",
     width = 4, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(0, 0, 0, 0), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(1, 2.8), ylim = c(-2, 2.2), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(1.7, -2, 1.7, 2)
segments(2.7, -2, 2.7, 2)
points(rep(1.7, 14), Tokyoso.ideal.point, pch = 19)
points(rep(2.7, 15), Zennorin.ideal.point, pch = 19)
segments(1.7, Tokyoso.ideal.point, 1.5, Tokyoso.ideal.point + Tokyoso.adjust)
segments(2.7, Zennorin.ideal.point, 2.5, Zennorin.ideal.point + Zennorin.adjust)
text(1.5, Tokyoso.ideal.point + Tokyoso.adjust, Tokyoso.labels, pos = 2)
text(2.5, Zennorin.ideal.point + Zennorin.adjust, Zennorin.labels, pos = 2)
points(1.73, median(Tokyoso.cutpoint), pch = 4)
text(1.73, median(Tokyoso.cutpoint), "cut-\npoint", pos = 4, cex = 0.8)
points(2.73, median(Zennorin.cutpoint), pch = 4)
text(2.73, median(Zennorin.cutpoint), "cut-\npoint", pos = 4, cex = 0.8)
text(1.35, 2.2, "Tokyo-to Kyoso case\n(April 2, 1969)", font = 2)
text(2.35, 2.2, "Zennorin Keishokuho case\n(April 25, 1973)", font = 2)
text(c(1.7, 2.7), -2, "\u2212", pos = 2)
text(c(1.7, 2.7), 2, "+", pos = 2)
dev.off()

#### Rationale for unidimensionality (analysis of median justices) ####
median.majority.sample <- matrix(NA, J, 6000)
other.majority.sample <- matrix(NA, J, 6000)
for (i in 1:J) {
  # Indicator for justices involved in the case
  participation <- (! is.na(data.matrix[i, ]))
  # Latent parameter draws of justices involved in the case in the year of the case
  ideal.point.sample.for.the.case <- ideal.point.sample.by.year[participation, as.numeric(nonunanimous.case.data$year[i]) - 1947, ]
  # Voting behavior in the case
  vote.in.the.case <- data.matrix[i, participation]
  # Conditional branching based on whether the number of justices is odd or even
  if (nonunanimous.case.data$n.justices[i] %% 2 == 1) {
    # Identify the median justice in each MCMC iteration
    who.is.the.median <- apply(ideal.point.sample.for.the.case, 2, function(x) which(x == median(x)))
    who.is.not.the.median <- t(apply(ideal.point.sample.for.the.case, 2, function(x) which(x != median(x))))
    # Record whether median and non-median justices took the majority side
    median.majority.sample[i, ] <- (vote.in.the.case[who.is.the.median] == 0) * 1
    for (j in 1:6000) {
      other.majority.sample[i, j] <- mean(vote.in.the.case[who.is.not.the.median[j, ]] == 0)
    }
  } else {
    # Identify the median justice in each MCMC iteration
    who.is.the.median <- t(apply(apply(ideal.point.sample.for.the.case, 2, rank), 2, 
                                 function(x) which(x == nrow(ideal.point.sample.for.the.case) / 2 | 
                                                     x == (nrow(ideal.point.sample.for.the.case) / 2) + 1)))
    who.is.not.the.median <- t(apply(apply(ideal.point.sample.for.the.case, 2, rank), 2, 
                                     function(x) which(x != nrow(ideal.point.sample.for.the.case) / 2 & 
                                                         x != (nrow(ideal.point.sample.for.the.case) / 2) + 1)))
    # Record whether median and non-median justices took the majority side
    for (j in 1:6000) {
      median.majority.sample[i, j] <- mean(vote.in.the.case[who.is.the.median[j, ]] == 0)
      other.majority.sample[i, j] <- mean(vote.in.the.case[who.is.not.the.median[j, ]] == 0)
    }
  }
}

# Probability of joining the majority for median justice(s)
round(mean(median.majority.sample), 3)
round(HPDinterval(mcmc(colMeans(median.majority.sample))), 3)

# Probability of joining the majority for non-median justices
round(mean(other.majority.sample), 3)
round(HPDinterval(mcmc(colMeans(other.majority.sample))), 3)

# Posterior probability that the former is smaller than the latter
round(mean(colMeans(median.majority.sample) - colMeans(other.majority.sample) < 0), 3)

#### Distribution of ideal points by career ####
## Posterior mean and credible interval of the average ideal point by career
career.labels <- c("Judge", "Attorney", "Prosecutor", "Bureaucrat", "Law Professor")
average.by.career <- matrix(NA, 5, 3)
for (i in 1:5) {
  average.by.career.sample <- rowMeans(ideal.point.whole[, justice.data$career == i])
  average.by.career[i, 1] <- mean(average.by.career.sample)
  average.by.career[i, 2:3] <- HPDinterval(mcmc(average.by.career.sample))
}
rownames(average.by.career) <- career.labels
round(average.by.career, 3)

# Calculate the posterior probability that the average for a career is smaller than that of another career (Footntoe 29)
career.comparison.matrix <- matrix(NA, 5, 5)
for (i in 1:5) {
  for (j in 1:5) {
    if (i == j) {
      next
    }
    career.comparison <- rowMeans(ideal.point.whole[, justice.data$career == i]) - 
      rowMeans(ideal.point.whole[, justice.data$career == j])
    career.comparison.matrix[i, j] <- sum(career.comparison < 0) / 6000
  }
}
round(career.comparison.matrix, 3)

## Distribution of ideal points for each career (Figure 3)
tiff("Figure_3.tiff",
     width = 4, height = 4.5, unit = "in", pointsize = 8, res = 1200, compression = "lzw")
layout(matrix(1:5, 5, 1))
par(mar = c(3, 11, 0, 0), lwd = 0.5)
for (i in c(3, 1, 4, 5, 2)) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(-4, 4), ylim = c(0, 1), 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = average.by.career[, 1], lty = 3, col = "gray80")
  hist(posterior.mean[justice.data$career == i], 
       freq = FALSE, col = "darkgray", border = NA, add = TRUE)
  segments(average.by.career[i, 2], 0.05, 
           average.by.career[i, 3], 0.05)
  points(average.by.career[i, 1], 0.05, pch = 19)
  axis(1, at = -4:4, lwd = 0.5)
  axis(2, at = seq(0, 0.8, 0.2), 
       labels = c("0.0", NA, "0.4", NA, "0.8"), lwd = 0.5)
  mtext(c("Liberal", "Conservative"), side = 1, line = 2, at = c(-4.05, 4.05), 
        adj = c(0, 1), cex = 0.6)
  mtext(paste0(career.labels[i], "\n(n = ", sum(justice.data$career == i), ")"), 
        line = -3, at = -6.8, adj = 0, 
        cex = 0.8, font = 2)
}
dev.off()

#### Analysis of justices who served as Chief Justice ####
## Comparison among all justices
# Count justices who served as Chief Justice
sum(justice.data$name %in% chief.justices)
# Calculate the average of the ideal points for justices who served as Chief Justice
round(mean(ideal.point.whole[, justice.data$name %in% chief.justices]), 3)
round(HPDinterval(mcmc(rowMeans(ideal.point.whole[, justice.data$name %in% chief.justices]))), 3)

# Count justices who did not serve as Chief Justice
sum(! justice.data$name %in% chief.justices)
# Calculate the average of the ideal points for justices who did not serve as Chief Justice
round(mean(ideal.point.whole[, ! justice.data$name %in% chief.justices]), 3)
round(HPDinterval(mcmc(rowMeans(ideal.point.whole[, ! justice.data$name %in% chief.justices]))), 3)

## Analysis of ideal points among career-judge justices
# Ideal points of career-judge justices
former.judge.ideal.point.whole <- ideal.point.whole[, justice.data$career == 1]
# Career-judge justices who served the Chief Justice
former.judge.chief.dummy <- justice.data$name[justice.data$career == 1] %in% chief.justices

# Count career-judge justices who served as Chief Justice
sum(former.judge.chief.dummy)
# Calculate the average of the ideal points for career-judge justices who served as Chief Justice
average.chief.est <- mean(former.judge.ideal.point.whole[, former.judge.chief.dummy])
average.chief.CI <- HPDinterval(mcmc(rowMeans(former.judge.ideal.point.whole[, former.judge.chief.dummy])))
round(average.chief.est, 3)
round(average.chief.CI, 3)

# Count career-judge justices who did not serve as Chief Justice
sum(! former.judge.chief.dummy)
# Calculate the average of the ideal points for career-judge justices who did not serve as Chief Justice
average.not.chief.est <- mean(former.judge.ideal.point.whole[, ! former.judge.chief.dummy])
average.not.chief.CI <- HPDinterval(mcmc(rowMeans(former.judge.ideal.point.whole[, ! former.judge.chief.dummy])))
round(average.not.chief.est, 3)
round(average.not.chief.CI, 3)

# Calculate the posterior probability that the average ideal points of career-judge justices who served as Chief Justice are greater than those who did not
former.judge.comparison <- rowMeans(former.judge.ideal.point.whole[, former.judge.chief.dummy]) - 
  rowMeans(former.judge.ideal.point.whole[, ! former.judge.chief.dummy])
round(sum(former.judge.comparison > 0) / 6000, 3)

## Plot the distribution of ideal points for career-judge justices, comparing those who served as Chief Justice to those who did not (Figure 4)
density.estimate.chief <- hist(posterior.mean[justice.data$career == 1 & 
                                                justice.data$name %in% chief.justices], 
                               plot = FALSE)
density.estimate.not.chief <- hist(posterior.mean[justice.data$career == 1 & 
                                                    ! justice.data$name %in% chief.justices], 
                                   plot = FALSE)

tiff("Figure_4.tiff",
     width = 4, height = 3, unit = "in", pointsize = 9, res = 1200, compression = "lzw")
par(mar = c(0, 1.5, 0, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-3, 3), ylim = c(-0.6, 1.2), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:length(density.estimate.chief$density)) {
  polygon(c(density.estimate.chief$breaks[i], density.estimate.chief$breaks[i], 
            density.estimate.chief$breaks[i + 1], density.estimate.chief$breaks[i + 1]), 
          c(0, density.estimate.chief$density[i], 
            density.estimate.chief$density[i], 0), 
          border = NA, col = "darkgray")
}
segments(average.chief.CI[1], 0.05, average.chief.CI[2], 0.05)
points(average.chief.est, 0.05, pch = 19)
for (i in 1:length(density.estimate.not.chief$density)) {
  polygon(c(density.estimate.not.chief$breaks[i], density.estimate.not.chief$breaks[i], 
            density.estimate.not.chief$breaks[i + 1], density.estimate.not.chief$breaks[i + 1]), 
          c(0, -1 * density.estimate.not.chief$density[i], 
            -1 * density.estimate.not.chief$density[i], 0), 
          border = NA, col = "darkgray")
}
segments(average.not.chief.CI[1], -0.05, average.not.chief.CI[2], -0.05)
points(average.not.chief.est, -0.05, pch = 19)
abline(h = 0)
segments(-3:3, 0.02, -3:3, -0.02)
text(-3:3, 0, -3:3, pos = 1, cex = 0.7)
axis(2, at = seq(-0.6, 1.2, 0.3), labels = NA, lwd = 0.5, cex.axis = 0.7)
mtext(c("0.6", "0.3", "0.0", "0.3", "0.6", "0.9", "1.2"), side = 2, 
      line = 0.75, at = seq(-0.6, 1.2, 0.3), cex = 0.7)
text(-3.2, 1.1, paste0("Who served the Chief Justice\n(n = ", 
                       sum(justice.data$career == 1 & 
                             justice.data$name %in% chief.justices), ")"), 
     pos = 4, font = 2)
text(-3.2, -0.5, paste0("Who did not serve\nthe Chief Justice\n(n = ", 
                        sum(justice.data$career == 1 & 
                              ! justice.data$name %in% chief.justices), ")"), 
     pos = 4, font = 2)
dev.off()

#### Analysis of the court median and voting patterns ####
## Modify justices' start and end years (see Footnote 30)
# Adjust the start and end years of each justice's tenure
justice.year.corrected <- cbind(as.numeric(format(justice.data$begin, "%Y")) - 1947, 
                                as.numeric(format(justice.data$end, "%Y")) - 1947)
justice.year.corrected[justice.year.corrected == 0] <- 1
justice.year.corrected[justice.year.corrected == 77] <- 76

ideal.point.sample.by.year.corrected <- ideal.point.sample.by.year
for (i in 1:N) {
  # Impute missing ideal points at the beginning of tenure for justices with incomplete data
  if (justice.year[i, 1] - justice.year.corrected[i, 1] == 1) {
    ideal.point.sample.by.year.corrected[i, justice.year[i, 1] - 1, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 1], ]
  } else if (justice.year[i, 1] - justice.year.corrected[i, 1] == 2) {
    ideal.point.sample.by.year.corrected[i, justice.year[i, 1] - 2, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 1] - 1, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 1], ]
  }
  # Exclude ideal points estimated for years after justices' retirement
  if (justice.year.corrected[i, 2] < 76) {
    ideal.point.sample.by.year.corrected[i, c((justice.year.corrected[i, 2] + 1):t), ] <- NA
  }
  # Impute missing ideal points at the end of tenure for justices with incomplete data
  if (justice.year[i, 2] - justice.year.corrected[i, 2] == -1) {
    ideal.point.sample.by.year.corrected[i, justice.year[i, 2] + 1, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 2], ]
  } else if (justice.year[i, 2] - justice.year.corrected[i, 2] == -2) {
    ideal.point.sample.by.year.corrected[i, justice.year[i, 2] + 2, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 2] + 1, ] <- 
      ideal.point.sample.by.year.corrected[i, justice.year[i, 2], ]
  }
}

## Track changes over time in the court's median ideal points (Figure 5)
# Identify active justices on December 31 for each year
Grand.Bench <- matrix(NA, N, t)
rownames(Grand.Bench) <- justice.data$english
for (i in 1:t) {
  Grand.Bench[, i] <- justice.data$begin < as.Date(paste0(1947 + i, "-12-31")) & 
    justice.data$end > as.Date(paste0(1947 + i, "-12-31"))
}

# Calculate the median ideal point of the court on December 31 annually
Grand.Bench.median <- matrix(NA, t, 3)
for (i in 1:t) {
  Grand.Bench.sample.by.year <- ideal.point.sample.by.year.corrected[Grand.Bench[, i], i, ]
  Grand.Bench.median.sample <- apply(Grand.Bench.sample.by.year, 2, median, na.rm = TRUE)
  Grand.Bench.median[i, 1] <- mean(Grand.Bench.median.sample, na.rm = TRUE)
  Grand.Bench.median[i, 2:3] <- HPDinterval(mcmc(Grand.Bench.median.sample))
}

tiff("Figure_5.tiff",
     width = 4, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(1, t), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(1, t, 10), t), col = "gray80")
abline(h = seq(-1, 1, 0.5), col = "gray80")
points(1:t, Grand.Bench.median[, 1], pch = 20)
segments(1:t, Grand.Bench.median[, 2], 1:t, Grand.Bench.median[, 3])
axis(1, at = c(seq(1, t, 10), t), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Year", 1, line = 2.5)
mtext("Median of ideal points", 2, line = 2.5)
dev.off()

## Analyze trends in discrimination parameters over time (Figure 6)
# Perform P-spline regression to model changes in discrimination parameters
gam.data <- data.frame(discrimination = discrimination.mean, 
                       date = as.numeric(nonunanimous.case.data$date), 
                       type = nonunanimous.case.data$type)
gam.result <- gam(discrimination ~ s(date, bs = "ps"), 
                  data = gam.data, method = "REML")
gam.prediction <- predict(gam.result, gam.data, type = "response")
gam.lpmatrix <- predict(gam.result, gam.data, type = "lpmatrix")
gam.se <- sqrt(rowSums((gam.lpmatrix %*% gam.result$Vc) * gam.lpmatrix))

tiff("Figure_6.tiff",
     width = 4, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", 
     xlim = range(as.numeric(nonunanimous.case.data$date)), ylim = c(-3, 3), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                "1978-01-01", "1988-01-01", "1998-01-01", 
                                "2008-01-01", "2018-01-01", "2023-01-01"))), 
       col = "gray80")
abline(h = -3:3, col = "gray80")
points(nonunanimous.case.data$date, discrimination.mean, 
       pch = ifelse(nonunanimous.case.data$n.justices > 5, 21, 23), col = NA, 
       bg = ifelse(discrimination.mean > 0, 
                   gray(apply(results.discrimination, 2, 
                              function(x) min(mean(x < 0) * 2, 1)), 0.5), 
                   gray(apply(results.discrimination, 2, 
                              function(x) min(mean(x > 0) * 2, 1)), 0.5)), 
       cex = ifelse(nonunanimous.case.data$n.justices > 5, 1, 0.8))
lines(gam.data$date, gam.prediction)
lines(gam.data$date, gam.prediction + qnorm(0.025) * gam.se, lty = 3)
lines(gam.data$date, gam.prediction + qnorm(0.975) * gam.se, lty = 3)
axis(1, at = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                  "1978-01-01", "1988-01-01", "1998-01-01", 
                                  "2008-01-01", "2018-01-01", "2023-01-01"))), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, at = -3:3, lwd = 0.5)
mtext("Date", 1, line = 2.5)
mtext("Discrimination parameter", 2, line = 2.5)
dev.off()

#### Analysis of appointer's partisanship ####
## Arrange the justices' data
# Flag justices appointed before the founding of the LDP; these data are not used in the analysis
justice.data$pre.LDP <- (justice.data$begin < "1955-11-15") * 1
# Justices appointed by non-LDP cabinets
justice.data$non.LDP <- ((justice.data$begin >= "1993-08-09" & 
                            justice.data$begin < "1996-01-11") | 
                           (justice.data$begin >= "2009-09-16" & 
                              justice.data$begin < "2012-12-26")) * 1
# Count the number of justices appointed after the LDP was founded
sum(justice.data$pre.LDP == 0)

## Estimation with Monte Carlo marginalization (see Appendix B.4)
cluster <- makeCluster(10)
registerDoParallel(cluster)
# Model (1)
lm.result.combined.1 <- foreach(i = 1:6000, .combine = "rbind", 
                                .packages = c("MCMCpack")) %dopar% {
                                  lm.data <- justice.data
                                  lm.data$theta <- ideal.point.whole[i, ]
                                  lm.result <- MCMCregress(theta ~ non.LDP + factor(career),
                                                           subset(lm.data, pre.LDP == 0),
                                                           burnin = 1000, mcmc = 1000,
                                                           seed = list(NA, i))
                                  lm.result
                                }
# Model (2)
lm.result.combined.2 <- foreach(i = 1:6000, .combine = "rbind", 
                                .packages = c("MCMCpack")) %dopar% {
                                  lm.data <- justice.data
                                  lm.data$theta <- ideal.point.whole[i, ]
                                  lm.result <- MCMCregress(theta ~ non.LDP + factor(career) +
                                                             non.LDP : I(career == 1),
                                                           subset(lm.data, pre.LDP == 0),
                                                           burnin = 1000, mcmc = 1000,
                                                           seed = list(NA, i))
                                  lm.result
                                }
# Model (3)
lm.result.combined.3 <- foreach(i = 1:6000, .combine = "rbind", 
                                .packages = c("MCMCpack")) %dopar% {
                                  lm.data <- justice.data
                                  lm.data$theta <- ideal.point.whole[i, ]
                                  lm.result <- MCMCregress(theta ~ non.LDP + factor(career) +
                                                             non.LDP : I(career != 2),
                                                           subset(lm.data, pre.LDP == 0),
                                                           burnin = 1000, mcmc = 1000,
                                                           seed = list(NA, i))
                                  lm.result
                                }
stopCluster(cluster)

## Results with Monte Carlo marginalization (Table 2)
# Model (1)
round(t(apply(lm.result.combined.1, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.combined.1, 2, posterior.summary)), 3)
round(mean(lm.result.combined.1[, 2] < 0), 3)
# Model (2)
round(t(apply(lm.result.combined.2, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.combined.2, 2, posterior.summary)), 3)
round(mean(lm.result.combined.2[, 2] < 0), 3)
# Model (3)
round(t(apply(lm.result.combined.3, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.combined.3, 2, posterior.summary)), 3)
round(mean(lm.result.combined.3[, 2] < 0), 3)
# 90% credible interval (Footnote 36)
round(HPDinterval(mcmc(lm.result.combined.3[, 2]), prob = 0.9), 3)

## Estimation without Monte Carlo marginalization (see Footnote 35 and Appendix C.5)
cluster <- makeCluster(10)
registerDoParallel(cluster)
# Model (1)
lm.result.alt.combined.1 <- foreach(i = 1:6000, .combine = "rbind", 
                                    .packages = c("MCMCpack")) %dopar% {
                                      lm.data <- justice.data
                                      lm.data$theta <- ideal.point.whole[i, ]
                                      lm.result <- lm(theta ~ non.LDP + factor(career),
                                                      lm.data, subset = pre.LDP == 0)
                                      lm.result$coefficients
                                    }
# Model (2)
lm.result.alt.combined.2 <- foreach(i = 1:6000, .combine = "rbind", 
                                    .packages = c("MCMCpack")) %dopar% {
                                      lm.data <- justice.data
                                      lm.data$theta <- ideal.point.whole[i, ]
                                      lm.result <- lm(theta ~ non.LDP + factor(career) +
                                                        non.LDP : I(career == 1),
                                                      lm.data, subset = pre.LDP == 0)
                                      lm.result$coefficients
                                    }
# Model (3)
lm.result.alt.combined.3 <- foreach(i = 1:6000, .combine = "rbind", 
                                    .packages = c("MCMCpack")) %dopar% {
                                      lm.data <- justice.data
                                      lm.data$theta <- ideal.point.whole[i, ]
                                      lm.result <- lm(theta ~ non.LDP + factor(career) +
                                                        non.LDP : I(career != 2),
                                                      lm.data, subset = pre.LDP == 0)
                                      lm.result$coefficients
                                    }
stopCluster(cluster)

## Results without Monte Carlo marginalization (Table A1)
# Model (1)
round(t(apply(lm.result.alt.combined.1, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.alt.combined.1, 2, posterior.summary)), 3)
round(mean(lm.result.alt.combined.1[, 2] < 0), 3)
# Model (2)
round(t(apply(lm.result.alt.combined.2, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.alt.combined.2, 2, posterior.summary)), 3)
round(mean(lm.result.alt.combined.2[, 2] < 0), 3)
# Model (3)
round(t(apply(lm.result.alt.combined.3, 2, posterior.summary, CI = FALSE)), 3)
round(t(apply(lm.result.alt.combined.3, 2, posterior.summary)), 3)
round(mean(lm.result.alt.combined.3[, 2] < 0), 3)

#### Comparison with previous studies (Appendix C.2)
## Hayakawa (1978)
# Extract latent trait parameters of justices served from 1955 to 1960 
ideal.point.sample.1955_60 <- ideal.point.sample.by.year[, 8:13, ]
target.justices.1955_60 <- ! is.na(rowMeans(ideal.point.sample.1955_60[ , , 1], na.rm = TRUE))
ideal.point.sample.1955_60 <- ideal.point.sample.1955_60[target.justices.1955_60, , ]
ideal.point.1955_60 <- apply(ideal.point.sample.1955_60, c(1, 3), mean, na.rm = TRUE)

# This study's estimates
posterior.mean.1955_60 <- rowMeans(ideal.point.1955_60)
names(posterior.mean.1955_60) <- justice.data$english[target.justices.1955_60]
posterior.CI.1955_60 <- HPDinterval(mcmc(t(ideal.point.1955_60)))

# Hayakawa's ranking
Hayakawa.rank <- justice.data$Hayakawa[target.justices.1955_60]

# Exclude justices with insufficient data from the analysis
target.justices.1955_60.limited <- 
  ! names(posterior.mean.1955_60) %in% c("Shimoyama, S", "Inoue, N", 
                                         "Takagi, T", "Ishisaka, S")
posterior.mean.1955_60 <- posterior.mean.1955_60[target.justices.1955_60.limited]
posterior.CI.1955_60 <- posterior.CI.1955_60[target.justices.1955_60.limited, ]
Hayakawa.rank <- rank(Hayakawa.rank[target.justices.1955_60.limited])

# Rank correlation coefficient
round(cor(posterior.mean.1955_60, Hayakawa.rank, method = "spearman"), 3)

# Figure A4
tiff("Figure_A4.tiff",
     width = 3, height = 2.5, unit = "in", pointsize = 6, res = 1200, compression = "lzw")
par(mar = c(3.5, 1, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-4, 4), ylim = c(1, 19), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-2, 1:19, 4, 1:19, lty = 3, col = "gray80")
segments(posterior.CI.1955_60[, 1], Hayakawa.rank, 
         posterior.CI.1955_60[, 2], Hayakawa.rank)
points(posterior.mean.1955_60, Hayakawa.rank, pch = 19)
text(-2.3, Hayakawa.rank, names(posterior.mean.1955_60), pos = 2)
axis(1, at = -2:4, lwd = 0.5)
mtext("This study's latent trait", side = 1, at = 1, line = 2.5)
mtext("Hayakawa's (1962) ranking", side = 2, line = 0)
dev.off()

## Osawa (1978)
# Extract latent trait parameters of justices served from 1960 to 1966 
ideal.point.sample.1960_66 <- ideal.point.sample.by.year[, 13:18, ]
target.justices.1960_66 <- ! is.na(rowMeans(ideal.point.sample.1960_66[ , , 1], na.rm = TRUE))
ideal.point.sample.1960_66 <- ideal.point.sample.1960_66[target.justices.1960_66, , ]
ideal.point.1960_66 <- apply(ideal.point.sample.1960_66, c(1, 3), mean, na.rm = TRUE)

# This study's estimates
posterior.mean.1960_66 <- rowMeans(ideal.point.1960_66)
names(posterior.mean.1960_66) <- justice.data$english[target.justices.1960_66]
posterior.CI.1960_66 <- HPDinterval(mcmc(t(ideal.point.1960_66)))

# Osawa's ranking
Osawa.rank <- justice.data$Osawa[target.justices.1960_66]

# Exclude justices with scarce data
target.justices.1960_66.limited <- 
  ! names(posterior.mean.1960_66) %in% c("Kotani, K", "Shima, T", "Saito, Y", 
                                         "Tanaka, K", "Saito, K", "Takahashi, K", 
                                         "Kusaka, A", "Osabe, K", "Ishida, K", 
                                         "Kido, Y", "Kashihara, G", "Tanaka, J", 
                                         "Matsuda, J", "Iwata, M", "Shimomura, S")
posterior.mean.1960_66 <- posterior.mean.1960_66[target.justices.1960_66.limited]
posterior.CI.1960_66 <- posterior.CI.1960_66[target.justices.1960_66.limited, ]
Osawa.rank <- rank(Osawa.rank[target.justices.1960_66.limited])

# Rank correlation coefficient
round(cor(posterior.mean.1960_66, Osawa.rank, method = "spearman"), 3)

# Figure A5
tiff("Figure_A5.tiff",
     width = 3, height = 2.5, unit = "in", pointsize = 6, res = 1200, compression = "lzw")
par(mar = c(3.5, 1, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-4.8, 2), ylim = c(1, 14), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-3, 1:14, 3, 1:14, lty = 3, col = "gray80")
segments(posterior.CI.1960_66[, 1], Osawa.rank, 
         posterior.CI.1960_66[, 2], Osawa.rank)
points(posterior.mean.1960_66, Osawa.rank, pch = 19)
text(-3.3, Osawa.rank, names(posterior.mean.1960_66), pos = 2)
axis(1, at = -3:2, lwd = 0.5)
mtext("This study's latent trait", side = 1, at = 0, line = 2.5)
mtext("Osawa's (1978) ranking", side = 2, line = 0)
dev.off()

## Comparison with Bushimata (1980)
# Extract latent trait parameters of justices served from 1967 to 1970
ideal.point.sample.1967_70 <- ideal.point.sample.by.year[, 19:22, ]
target.justices.1967_70 <- ! is.na(rowMeans(ideal.point.sample.1967_70[ , , 1], na.rm = TRUE))
ideal.point.sample.1967_70 <- ideal.point.sample.1967_70[target.justices.1967_70, , ]
ideal.point.1967_70 <- apply(ideal.point.sample.1967_70, c(1, 3), mean, na.rm = TRUE)

# This study's estimates
posterior.mean.1967_70 <- rowMeans(ideal.point.1967_70)
names(posterior.mean.1967_70) <- justice.data$english[target.justices.1967_70]
posterior.CI.1967_70 <- HPDinterval(mcmc(t(ideal.point.1967_70)))

# Bushimata's score
Bushimata.score <- justice.data$Bushimata[target.justices.1967_70] * -1

# Correlation coefficient
round(cor(posterior.mean.1967_70, Bushimata.score, use = "p"), 3)

# Figure A6
tiff("Figure_A6.tiff",
     width = 3, height = 2.5, unit = "in", pointsize = 6, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-2, 2), ylim = c(-0.5, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(posterior.CI.1967_70[, 1], Bushimata.score, 
         posterior.CI.1967_70[, 2], Bushimata.score)
points(posterior.mean.1967_70, Bushimata.score, pch = 19)
text(posterior.mean.1967_70, Bushimata.score - 0.02, 
     names(posterior.mean.1967_70), pos = 3, cex = 0.7)
axis(1, at = -2:2, lwd = 0.5)
axis(2, at = seq(-0.5, 1, 0.5), lwd = 0.5)
mtext("This study's latent trait", side = 1, line = 2.5)
mtext("Bushimata's (1980) score (inverted)", side = 2, line = 2.5)
dev.off()

## Itoh (2010)
# extract latent trait parameters of justices served from 1948 to 2007
ideal.point.sample.1948_2007 <- ideal.point.sample.by.year[, 1:60, ]
target.justices.1948_2007 <- ! is.na(rowMeans(ideal.point.sample.1948_2007[ , , 1], na.rm = TRUE))
ideal.point.sample.1948_2007 <- ideal.point.sample.1948_2007[target.justices.1948_2007, , ]
ideal.point.1948_2007 <- apply(ideal.point.sample.1948_2007, c(1, 3), mean, na.rm = TRUE)

# This study's estimates
posterior.mean.1948_2007 <- rowMeans(ideal.point.1948_2007)
names(posterior.mean.1948_2007) <- justice.data$english[target.justices.1948_2007]
posterior.CI.1948_2007 <- HPDinterval(mcmc(t(ideal.point.1948_2007)))

# Itoh's score
justice.data$Itoh <- justice.data$con * 2 / 
  (justice.data$lib + justice.data$con)
# Exclude justices with scarce data and with an error in Itoh's list
justice.data$Itoh[justice.data$lib + justice.data$con < 5] <- NA
justice.data$Itoh[justice.data$english == "Kobayashi, S"] <- NA
Itoh.score <- justice.data$Itoh[target.justices.1948_2007]
Itoh.score[justice.data$lib[target.justices.1948_2007] + 
             justice.data$con[target.justices.1948_2007] < 5] <- NA

# Number of observations
sum(! is.na(Itoh.score))

# Correlation coefficient
round(cor(posterior.mean.1948_2007, Itoh.score, use = "p"), 3)

# Figure A7
tiff("Figure_A7.tiff",
     width = 3, height = 2.5, unit = "in", pointsize = 6, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-3, 4), ylim = c(0, 2), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(posterior.CI.1948_2007[, 1], Itoh.score, 
         posterior.CI.1948_2007[, 2], Itoh.score)
points(posterior.mean.1948_2007, Itoh.score, pch = 19)
axis(1, at = -3:4, lwd = 0.5)
axis(2, at = seq(0, 2, 0.5), lwd = 0.5)
mtext("This study's latent trait", side = 1, line = 2.5)
mtext("Itoh's (2010) score (inverted)", side = 2, line = 2.5)
dev.off()

#### Analysis of the court median based on Itoh's score (Appendix C.3) ####
# Court median based on Itoh's score
Grand.Bench.median.Itoh <- rep(NA, t)
for (i in 1:t) {
  Grand.Bench.justices.by.year <- justice.data$Itoh[Grand.Bench[, i]]
  if (sum(is.na(Grand.Bench.justices.by.year)) < 8) {
    Grand.Bench.median.Itoh[i] <- median(Grand.Bench.justices.by.year, na.rm = TRUE)
  }
}

# Temporal change in the court median based on Itoh's score (Figure A8)
tiff("Figure_A8.tiff",
     width = 4, height = 3, unit = "in", pointsize = 7, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 3.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(1, t), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(1, t, 10), t), col = "gray80")
abline(h = seq(-1, 1, 0.5), col = "gray80")
points(1:t, Grand.Bench.median[, 1], pch = 20, col = "darkgray")
segments(1:t, Grand.Bench.median[, 2], 1:t, Grand.Bench.median[, 3], col = "darkgray")
lines(1:t, Grand.Bench.median.Itoh - 1, lwd = 1)
axis(1, at = c(seq(1, t, 10), t), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, at = seq(-1, 1, 0.5), labels = seq(0, 2, 0.5), lwd = 0.5)
axis(4, lwd = 0.5)
mtext("Year", 1, line = 2.5)
mtext("Median of Itoh's (2010) score (solid line)", 2, line = 2.5)
mtext("Median of ideal points (dots and segments)", 4, line = 2.5)
dev.off()

#### Analysis of voting patterns by case category (Appendix C.4) ####
# P-spline regression for civil cases
gam.result.civil <- gam(discrimination ~ s(date, bs = "ps"), data = gam.data, 
                        subset = type == "civil", method = "REML")
gam.prediction.civil <- predict(gam.result.civil, 
                                subset(gam.data, type == "civil"), 
                                type = "response")
gam.lpmatrix.civil <- predict(gam.result.civil, 
                              subset(gam.data, type == "civil"), 
                              type = "lpmatrix")
gam.se.civil <- sqrt(rowSums((gam.lpmatrix.civil %*% gam.result.civil$Vc) * 
                               gam.lpmatrix.civil))

# P-spline regression for criminal cases
gam.result.criminal <- gam(discrimination ~ s(date, bs = "ps"), data = gam.data, 
                           subset = type == "criminal", method = "REML")
gam.prediction.criminal <- predict(gam.result.criminal, 
                                   subset(gam.data, type == "criminal"), 
                                   type = "response")
gam.lpmatrix.criminal <- predict(gam.result.criminal, 
                                 subset(gam.data, type == "criminal"), 
                                 type = "lpmatrix")
gam.se.criminal <- sqrt(rowSums((gam.lpmatrix.criminal %*% gam.result.criminal$Vc) * 
                                  gam.lpmatrix.criminal))

# P-spline regression for administrative cases
gam.result.administrative <- gam(discrimination ~ s(date, bs = "ps"), data = gam.data, 
                                 subset = type == "administrative", method = "REML")
gam.prediction.administrative <- predict(gam.result.administrative, 
                                         subset(gam.data, type == "administrative"), 
                                         type = "response")
gam.lpmatrix.administrative <- predict(gam.result.administrative, 
                                       subset(gam.data, type == "administrative"), 
                                       type = "lpmatrix")
gam.se.administrative <- sqrt(rowSums((gam.lpmatrix.administrative %*% gam.result.administrative$Vc) * 
                                        gam.lpmatrix.administrative))

# Trends in discrimination parameters over time by case category (Figure A9)
tiff("Figure_A9.tiff",
     width = 6, height = 4.5, unit = "in", pointsize = 8, res = 1200, compression = "lzw")
layout(matrix(1:4, 2, 2, byrow = TRUE))
par(mar = c(3.5, 3.5, 2, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = range(as.numeric(nonunanimous.case.data$date)), ylim = c(-3, 3), 
     main = "Civil cases", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                "1978-01-01", "1988-01-01", "1998-01-01", 
                                "2008-01-01", "2018-01-01", "2023-01-01"))), 
       col = "gray80")
abline(h = -3:3, col = "gray80")
points(nonunanimous.case.data$date[nonunanimous.case.data$type == "civil"], 
       discrimination.mean[nonunanimous.case.data$type == "civil"], 
       pch = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "civil"] > 5, 21, 23), col = NA, 
       bg = ifelse(discrimination.mean[nonunanimous.case.data$type == "civil"] > 0, 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "civil"], 2, 
                              function(x) min(mean(x < 0) * 2, 1)), 0.5), 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "civil"], 2, 
                              function(x) min(mean(x > 0) * 2, 1)), 0.5)), 
       cex = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "civil"] > 5, 1, 0.8))
lines(gam.data$date[gam.data$type == "civil"], 
      gam.prediction.civil)
lines(gam.data$date[gam.data$type == "civil"], 
      gam.prediction.civil + qnorm(0.025) * gam.se.civil, lty = 3)
lines(gam.data$date[gam.data$type == "civil"], 
      gam.prediction.civil + qnorm(0.975) * gam.se.civil, lty = 3)
axis(1, at = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                  "1978-01-01", "1988-01-01", "1998-01-01", 
                                  "2008-01-01", "2018-01-01", "2023-01-01"))), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, at = -3:3, lwd = 0.5)
mtext("Date", 1, line = 2.5)
mtext("Discrimination parameter", 2, line = 2.5)
plot(NULL, NULL, bty = "n", xlim = range(as.numeric(nonunanimous.case.data$date)), ylim = c(-3, 3), 
     main = "Criminal cases", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                "1978-01-01", "1988-01-01", "1998-01-01", 
                                "2008-01-01", "2018-01-01", "2023-01-01"))), 
       col = "gray80")
abline(h = -3:3, col = "gray80")
points(nonunanimous.case.data$date[nonunanimous.case.data$type == "criminal"], 
       discrimination.mean[nonunanimous.case.data$type == "criminal"], 
       pch = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "criminal"] > 5, 21, 23), col = NA, 
       bg = ifelse(discrimination.mean[nonunanimous.case.data$type == "criminal"] > 0, 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "criminal"], 2, 
                              function(x) min(mean(x < 0) * 2, 1)), 0.5), 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "criminal"], 2, 
                              function(x) min(mean(x > 0) * 2, 1)), 0.5)), 
       cex = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "criminal"] > 5, 1, 0.8))
lines(gam.data$date[gam.data$type == "criminal"], 
      gam.prediction.criminal)
lines(gam.data$date[gam.data$type == "criminal"], 
      gam.prediction.criminal + qnorm(0.025) * gam.se.criminal, lty = 3)
lines(gam.data$date[gam.data$type == "criminal"], 
      gam.prediction.criminal + qnorm(0.975) * gam.se.criminal, lty = 3)
axis(1, at = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                  "1978-01-01", "1988-01-01", "1998-01-01", 
                                  "2008-01-01", "2018-01-01", "2023-01-01"))), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, at = -3:3, lwd = 0.5)
mtext("Date", 1, line = 2.5)
mtext("Discrimination parameter", 2, line = 2.5)
plot(NULL, NULL, bty = "n", xlim = range(as.numeric(nonunanimous.case.data$date)), ylim = c(-3, 3), 
     main = "Administrative cases", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                "1978-01-01", "1988-01-01", "1998-01-01", 
                                "2008-01-01", "2018-01-01", "2023-01-01"))), 
       col = "gray80")
abline(h = -3:3, col = "gray80")
points(nonunanimous.case.data$date[nonunanimous.case.data$type == "administrative"], 
       discrimination.mean[nonunanimous.case.data$type == "administrative"], 
       pch = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "administrative"] > 5, 21, 23), col = NA, 
       bg = ifelse(discrimination.mean[nonunanimous.case.data$type == "administrative"] > 0, 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "administrative"], 2, 
                              function(x) min(mean(x < 0) * 2, 1)), 0.5), 
                   gray(apply(results.discrimination[, nonunanimous.case.data$type == "administrative"], 2, 
                              function(x) min(mean(x > 0) * 2, 1)), 0.5)), 
       cex = ifelse(nonunanimous.case.data$n.justices[nonunanimous.case.data$type == "administrative"] > 5, 1, 0.8))
lines(gam.data$date[gam.data$type == "administrative"], 
      gam.prediction.administrative)
lines(gam.data$date[gam.data$type == "administrative"], 
      gam.prediction.administrative + qnorm(0.025) * gam.se.administrative, lty = 3)
lines(gam.data$date[gam.data$type == "administrative"], 
      gam.prediction.administrative + qnorm(0.975) * gam.se.administrative, lty = 3)
axis(1, at = as.numeric(as.Date(c("1948-01-01", "1958-01-01", "1968-01-01", 
                                  "1978-01-01", "1988-01-01", "1998-01-01", 
                                  "2008-01-01", "2018-01-01", "2023-01-01"))), 
     labels = c(seq(48, 98, 10), "08", "18", "23"), lwd = 0.5)
axis(2, at = -3:3, lwd = 0.5)
mtext("Date", 1, line = 2.5)
mtext("Discrimination parameter", 2, line = 2.5)
dev.off()

#### Robustness check excluding cases coded solely by the author (Appendix C.6) ####
# Count the number of cases that were coded solely by the author
table(voting.records$author)
# Count the cases where the author's original codings were maintained or new codings different from other coders were adopted
table(voting.records$filename[voting.records$author == 1])

# Exclude cases that were coded solely by the author to prepare the data matrix for analysis
data.matrix.r <- data.matrix[voting.records$author == 0, ]
nonunanimous.case.data.r <- nonunanimous.case.data[voting.records$author == 0, ]

## Preliminary estimation of the static IRT model
cluster <- makeCluster(10)
registerDoParallel(cluster)
static.results.r <- foreach(i = 1:10, .packages = "MCMCpack") %dopar% {
  data.matrix.sub <- data.matrix.r[nonunanimous.case.data.r$year >= start.year[i] &
                                     nonunanimous.case.data.r$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  data.matrix.sub <- data.matrix.sub[, target.justices.sub]
  MCMCirt1d(datamatrix = t(data.matrix.sub), 
            burnin = 100000, mcmc = 100000, thin = 50, 
            seed = list(NA, i), AB0 = matrix(c(0.64, 0, 0, 0.64), 2, 2))
}
stopCluster(cluster)

# save(static.results.r, file = "static_results_r.Rdata")
# load("static_results_r.Rdata")

## Estimated latent trait parameters for each period
layout(matrix(1:10, 5, 2, byrow = TRUE))
par(mar = c(2.5, 0.5, 2, 0.5), lwd = 0.25)
for (i in 1:10) {
  data.matrix.sub <- data.matrix.r[nonunanimous.case.data.r$year >= start.year[i] &
                                     nonunanimous.case.data.r$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  justice.data.sub <- justice.data[target.justices.sub, ]
  results.theta <- static.results[[i]][, stri_detect_regex(colnames(static.results.r[[i]]), "theta.")]
  iteration.mean <- rowMeans(results.theta)
  iteration.sd <- apply(results.theta, 1, sd)
  for (j in 1:nrow(results.theta)) {
    results.theta[j, ] <- (results.theta[j, ] - iteration.mean[j]) / iteration.sd[j]
  }
  posterior.mean <- colMeans(results.theta)
  posterior.CI <- HPDinterval(mcmc(results.theta), 0.95)
  posterior.order <- order(posterior.mean)
  justice.data.sub$pch <- ifelse(justice.data.sub$career == 1, 21, 
                                 ifelse(justice.data.sub$career == 2, 22, 
                                        ifelse(justice.data.sub$career == 3, 4, 
                                               ifelse(justice.data.sub$career == 4, 24, 23))))
  justice.data.sub$bg <- c("black", "white", NA, "darkgray", "lightgray")[justice.data.sub$career]
  plot(NULL, NULL, bty = "n", xlim = c(-4.5, 4.5), ylim = c(1, length(posterior.mean)),
       main = paste0(start.year[i], "\u2013", end.year[i]),
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = -4:4, col = "gray")
  segments(posterior.CI[posterior.order, 1], c(1:nrow(posterior.CI)),
           posterior.CI[posterior.order, 2], c(1:nrow(posterior.CI)))
  points(posterior.mean[posterior.order], 1:length(posterior.mean),
         pch = justice.data.sub$pch[posterior.order],
         bg = justice.data.sub$bg[posterior.order], cex = 0.6)
  text(ifelse(rep(c(0, 1), 100)[1:length(posterior.mean)] == 0,
              posterior.CI[posterior.order, 1],
              posterior.CI[posterior.order, 2]), 1:length(posterior.mean),
       justice.data.sub$english[posterior.order],
       pos = rep(c(2, 4), 100)[1:length(posterior.mean)], cex = 0.6,
       font = ifelse(justice.data.sub$name[posterior.order] %in% chief.justices, 2, 1))
  axis(1, at = -4:4, lwd = 0.25)
}

## Identify and list the most extreme justices based on their latent trait estimates for each period
extreme.judges.r.jp <- extreme.judges.r.en <- list()
for (i in 1:10) {
  data.matrix.sub <- data.matrix.r[nonunanimous.case.data.r$year >= start.year[i] &
                                     nonunanimous.case.data.r$year <= end.year[i], ]
  target.justices.sub <- apply(data.matrix.sub, 2, function(x) sum(! is.na(x))) > 0
  justice.data.sub <- justice.data[target.justices.sub, ]
  results.theta <- static.results.r[[i]][, stri_detect_regex(colnames(static.results.r[[i]]), "theta.")]
  posterior.mean <- colMeans(results.theta)
  posterior.order <- order(posterior.mean)
  extreme.judges.r.jp[[i]] <- as.character(justice.data.sub$name[posterior.order][c(1, length(posterior.order))])
  extreme.judges.r.en[[i]] <- as.character(justice.data.sub$english[posterior.order][c(1, length(posterior.order))])
}
extreme.judges.r.en

# Signs of the latent traits are fixed in the estimation of the dynamic IRT model
# The above figure tells us which side in a period is connected to which side in the next period.
sign.r <- c(-1, 1, 1, 1, -1, -1, -1, 1, -1, -1)
theta.constraints.r <- list()
for (i in 1:10) {
  if (sign.r[i] == 1) {
    theta.constraints.r[[2 * i - 1]] <- "-"
    theta.constraints.r[[2 * i]] <- "+"
  } else {
    theta.constraints.r[[2 * i - 1]] <- "+"
    theta.constraints.r[[2 * i]] <- "-"
  }
  names(theta.constraints.r)[2 * i - 1] <- extreme.judges.r.jp[[i]][1]
  names(theta.constraints.r)[2 * i] <- extreme.judges.r.jp[[i]][2]
}

## Estimation of the dynamic IRT model
# Year id
year.r <- as.numeric(nonunanimous.case.data.r$year) - 
  as.numeric(min(nonunanimous.case.data.r$year)) + 1
# Number of cases
J.r <- nrow(data.matrix.r)

cluster <- makeCluster(10)
registerDoParallel(cluster)
dynamic.list.r <- foreach(i = 1:3, .packages = "MCMCpack") %dopar% {
  MCMCdynamicIRT1d(datamatrix = t(data.matrix.r), 
                   item.time.map = year.r, 
                   theta.constraints = theta.constraints.r, 
                   burnin = 200000, mcmc = 200000, 
                   thin = 100, seed = list(NA, i), 
                   tau2.start = 0.1, A0 = 0.64, B0 = 0.64)
}
stopCluster(cluster)

# save(dynamic.list.r, file = "dynamic_list_r.Rdata")
# load("dynamic_list_r.Rdata")

# Convergence check
dynamic.results.r <- mcmc.list(dynamic.list.r)
which(gelman.diag(dynamic.results.r, multivariate = FALSE)$psrf[, 1] > 1.1)

## Post-processing of the posterior draws
raw.MCMC.draws.r <- normalized.MCMC.draws.r <- 
  rbind(dynamic.results.r[[1]], dynamic.results.r[[2]], dynamic.results.r[[3]])
iteration.mean.r <- rowMeans(raw.MCMC.draws[, initial.period.column])
iteration.sd.r <- apply(raw.MCMC.draws[, initial.period.column], 1, sd)
for (i in 1:nrow(raw.MCMC.draws)) {
  normalized.MCMC.draws.r[i, stri_detect_regex(colnames(normalized.MCMC.draws.r), "theta")] <- 
    (raw.MCMC.draws.r[i, stri_detect_regex(colnames(raw.MCMC.draws.r), "theta")] - iteration.mean.r[i]) / iteration.sd.r[i]
  normalized.MCMC.draws.r[i, stri_detect_regex(colnames(normalized.MCMC.draws.r), "alpha")] <- 
    raw.MCMC.draws.r[i, stri_detect_regex(colnames(raw.MCMC.draws.r), "alpha")] - 
    iteration.mean.r[i] * raw.MCMC.draws.r[i , stri_detect_regex(colnames(raw.MCMC.draws.r), "beta")]
  normalized.MCMC.draws.r[i , stri_detect_regex(colnames(normalized.MCMC.draws.r), "beta")] <- 
    iteration.sd.r[i] * raw.MCMC.draws.r[i , stri_detect_regex(colnames(raw.MCMC.draws.r), "beta")]
}
normalized.MCMC.draws.r <- mcmc(normalized.MCMC.draws.r)

## Posterior draws for each justice's average over their tenure
ideal.point.whole.r <- matrix(NA, nrow(normalized.MCMC.draws.r), N)
for (i in 1:N) {
  temp.draws <- normalized.MCMC.draws.r[, stri_detect_regex(colnames(normalized.MCMC.draws.r), paste0("theta.", justice.data$name[i]))]
  if (is.null(dim(temp.draws))) {
    ideal.point.whole.r[, i] <- temp.draws
  } else {
    ideal.point.whole.r[, i] <- rowMeans(temp.draws)
  }
}

## Comparison with the original results (Figure A10)
tiff("Figure_A10.tiff",
     width = 4, height = 4, unit = "in", pointsize = 8, res = 1200, compression = "lzw")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 3.5), ylim = c(-2.5, 3.5), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = -2:3, col = "lightgray")
abline(h = -2:3, col = "lightgray")
box()
points(colMeans(ideal.point.whole), colMeans(ideal.point.whole.r))
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Original estimates", side = 1, line = 2.5)
mtext("Estimates using data excluding cases decided by the author's judgment", 
      side = 2, line = 2.5)
dev.off()

# Correlation with the original results
round(cor(colMeans(ideal.point.whole), colMeans(ideal.point.whole.r)), 3)